mapSNPsToContigs <- function(contig.dataframe,contig.map,cDNA.name){
	#this function extracts snps from individual sequence captured contigs
	#and maps them onto layout corresponding to the structure of 
	#the join.contigs function
	#
	#contig.dataframe is produced by readNewblerMapper
	#contig.map is produced by join.contigs
	#
	tmp.snp.out <- list()
	for(i in 1:length(contig.map[[cDNA.name]]$sc.contig)){
		gctg <- contig.map[[cDNA.name]]$sc.contig[i]
	
		tmp.snp.in <- grouse1hc[grep(gctg,contig.dataframe$ref.accno),]
	
		#jump to next i if no snps
		if(dim(tmp.snp.in)[1]==0) next
	
		tmp.concat.end <- c(0,contig.map[[cDNA.name]]$concat.end)
		tmp.snp.out[[gctg]] <- tmp.snp.in
	
		#change start and ends relative to matched contig, and reverse comp if needed
		if(contig.map[[cDNA.name]]$direction[i]=="forward"){
			tmp.snp.out[[gctg]]$start.pos <- tmp.concat.end[i] + tmp.snp.in$start.pos
			tmp.snp.out[[gctg]]$end.pos <- tmp.concat.end[i] + tmp.snp.in$end.pos
		
			}
		if(contig.map[[cDNA.name]]$direction[i]=="reverse"){
			cat('reversing \n')
			tmp.snp.out[[gctg]]$start.pos <- tmp.concat.end[i] + contig.map[[cDNA.name]]$length[i] - tmp.snp.in$start.pos + 1
		tmp.snp.out[[gctg]]$end.pos <- tmp.concat.end[i] + contig.map[[cDNA.name]]$length[i] - tmp.snp.in$end.pos + 1
			tmp.snp.out[[gctg]]$ref.nucl[tmp.snp.in$ref.nucl!="-"] <- 
				sapply(tmp.snp.in$ref.nucl[tmp.snp.in$ref.nucl!="-"],strComp)
			tmp.snp.out[[gctg]]$var.nucl[tmp.snp.in$var.nucl!="-"] <- 
				sapply(tmp.snp.in$var.nucl[tmp.snp.in$var.nucl!="-"],strComp)

			}
		if(contig.map[[cDNA.name]]$direction[i] %in% c("forward","reverse") == FALSE){
			warning('neither forward nor reverse for', gctg, "in", names(contig.map)[cDNA.name])
			}
		#output is a list of dataframes, one per sequence captured contig
		}
 
 	#tidy up list into a single dataframe
	 if(length(tmp.snp.out)>0){
 		tmp.snp.df <- tmp.snp.out[[1]]
 		if(length(tmp.snp.out)>1){
 			for(i in 2:length(tmp.snp.out)){
 				tmp.snp.df <- rbind(tmp.snp.df,tmp.snp.out[[i]])
 				}
 		
 			}
 		}
	tmp.snp.df
	}

#can't get the function below to work
mapCoverageQual <- function(AceObject,qual.location,qual.threshold){
	#takes input from makeAceObject
	#and location of qual files
	#to calculate coverage for each base
	#above qual.threshold
	#output is a list
	
	if(qual.threshold>40 || qual.threshold<0) warning('quality threshold not in [0,40]')
	
	tmp.pad.df <- data.frame(pad.pos=1:nchar(AceObject$padded.seq),
		base=strsplit(AceObject$padded.seq,"")[[1]],
		coverage=0,stringsAsFactors=FALSE)
	
	tmp.AF <- AceObject$AF[-grep('contig',AceObject$AF$read),]
	tmp.QA <- AceObject$QA[-grep('contig',AceObject$QA$read),]
	#tmp.RD <- AceObject$QA[-grep('contig',AceObject$RD$read),]
	tmp.read.seq <- AceObject$read.seq[-grep('contig',names(AceObject$read.seq))]
	
	#tmp.read.len <- tmp.QA$stop.pos.clear + 1 - tmp.QA$start.pos.clear
	tmp.cont.stop <- tmp.AF$start.pos + tmp.QA$stop.pos.clear -1
	tmp.cont.start <- tmp.AF$start.pos + tmp.QA$start.pos.clear -1
	tmp.cont.start[tmp.cont.start < 1] <- 1

	tmp.filter <- tmp.cont.stop > 1
	tmp.cont.stop <- tmp.cont.stop[tmp.filter]
	tmp.cont.start <- tmp.cont.start[tmp.filter]

	if(!identical(tmp.AF$read,tmp.QA$read)) error('AF does not match RD')
	
	#read in quality info
	tmp.qf <- readQualFiles(qual.reads=tmp.QA$read,location=qual.location)
	cat('read quality info for',paste(tmp.QA$read,collapse=", "),'\n')
	tmp.qual <- extractQual(tmp.qf)
	
	qual.padded <- list()
	for(i in 1:length(tmp.qual)){
		tmp.qrd <- tmp.qual[[i]]
		if(tmp.AF$orientation[grep(names(tmp.qual)[i],tmp.AF$read)[1]]=="C"){ #complement
			tmp.qrd <- tmp.qrd[length(tmp.qrd):1]
			}
		tmp.rdi <- strsplit(tmp.read.seq[[i]],"")[[1]]
		qual.padded[[i]] <- rep(0,length(tmp.rdi))  
		qual.padded[[i]][-grep("\\*",tmp.rdi)] <- tmp.qrd
		}
	names(qual.padded) <- names(tmp.read.seq)
	for(i in 1:length(tmp.cont.stop)){
		#add 1 to coverage for each base for each overlapping read, weighted by quality
		
		#qual.vec <- rep(0,nchar(AceObject$padded.seq))
		#qual.map <- tmp.pad.df$pad.pos[-grep("\\*",tmp.pad.df$base)]
		
		
		tmp.pad.df$coverage[tmp.cont.start[i]:tmp.cont.stop[i]] <- tmp.pad.df$coverage[tmp.cont.start[i]:tmp.cont.stop[i]] +1
		}
	list(coverage=tmp.pad.df[tmp.pad.df$base!="*",],padded.coverage=tmp.pad.df)
	}
 
 